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Abstract. Preferential attachment is a popular model of growing networks. We consider a generalized 
model with random node removal, and a combination of preferential and random attachment. Using a 
high-degree expansion of the master equation, we identify a topological phase transition depending on the 
rate of node removal and the relative strength of preferential vs. random attachment, where the degree 
distribution goes from a power law to one with an exponential tail. 



1 Introduction 

Complex networks are found in nature, social and eco- 
nomic systems, technical infrastructures, and countless 
other fields. The macroscopic properties of such networks 
emerge from the microscopic interaction of many individ- 
ual constituents. Various models of complex networks have 
been proposed and the statistical mechanics of networks 
has become an established branch of statistical physics 
[1,2,3,4,5,6]. Since complex networks are non-equilibrium 
systems, they do not have to obey detailed balance, and 
may show many fascinating features not found in equilib- 
rium systems. 

In order to explain the power-law degree distribution 
observed in many complex networks, Barabasi and Al- 
bert [7] introduced a preferential attachment model for 
growing networks. When new nodes enter the network, 
they prefer to attach to nodes with high degree. In a gener- 
alization of this model, the probability p a dd(v) that a new 
node u establishes an edge to an existing node v with de- 
gree k(v) is proportional to an attractiveness function A^; 
normalizing, _p a dd(«) = ^-k(v)/Y, w A k(w)- Barabasi and 
Albert considered the case where the attachment is pro- 
portional to degree, Ak = k. This results in a power-law 
degree distribution ~ k f with exponent 7 = 3, which 
is close to the observed exponents of many real world net- 
works [3]. 

However, the linear preferential attachment function 
of the Barabasi-Albert model was introduced as an ad 
hoc ansatz without fundamental justification, and many 
generalizations are conceivable. For many networks, a re- 



alistic model has to take into account node removal, edge 
rewiring [8] or removal [9], and other dynamical processes, 
as well as deviations from linear preferential attachment. 

Here we consider a generalized preferential attachment 
model with asymptotically linear attractiveness function 
and random node removal. 

The degree distribution of generalized preferential at- 
tachment models is very sensitive to model-specific fea- 
tures. Varying the attractiveness function [10,11,12,13,14] 
or including node removal [15,16,17] can shift the expo- 
nent of the power-law degree distribution to 2 < 7 < 00. 
For networks of constant size [17] or a sublinear attrac- 
tiveness function [11], the degree distribution can become 
a stretched exponential. We will demonstrate that gener- 
alizations of the Barabasi-Albert model can dramatically 
affect the degree distribution even in the case of growing 
networks and asymptotically linear attractiveness, leading 
to a topological phase transition from a power-law degree 
distribution to an exponential degree distribution, with a 
stretched exponential at the critical point. 

2 The Model 

The topological phase transition in generalized preferen- 
tial attachment networks can be illustrated by consider- 
ing the following model. Vertices arrive at rate 1, each 
new vertex makes c connections to existing vertices, and 
we remove vertices randomly at a rate r. Each new edge 
attaches to a given pre-existing vertex of degree k with 
probability proportional to its attractiveness A k . We as- 
sume that Ak is of the form 
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A k = k + k 
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for some constant k* . Thus the choice of a link endpoint 
is somewhere between preferential attachment (the case 
k* = 0) and uniform attachment (the limit k* — > oo). We 
can also treat k* as a the initial degree of the vertex when 
it is added to the network (e.g. by adding k* self-loops) 
and then run the "pure" preferential attachment model. 

The network starts to grow at time t = with No 
nodes and Mq links. We are interested in systems where 
the number of nodes is much larger than unity. Since we 
want to study the influence on the topology of the dy- 
namics of our model, and not the influence of the initial 
network, we let the system evolve until we have added a 
number of nodes much larger than N , and the degree dis- 
tribution has reached equilibrium. For r < 1 the number 
of nodes grows with time, so the initial value Nq is unim- 
portant. For r = 1, nodes are added and removed at the 
same rate; in that case the expected number of nodes is 
constant, so we start with No 3> 1. 



3 Analytic Solution 

In this section, we derive the average degree (k) and the 
average attractiveness (A) . We then write the master equa- 
tion for the degree distribution pk and solve for its asymp- 
totic behavior for large k in terms of the parameters r, c, 
and k*. 



3.1 Mean Degree and Attractiveness 

Let p k be the expected fraction of vertices in the network 
at a given time that have degree k. As in [17], the expected 
mean degree of a vertex (k) = J2T=o ^Pk can ^ e derived as 
follows. The expected increase in the number of vertices 
per unit time is 1 — r. The expected number of edges re- 
moved when a randomly chosen vertex is removed is (fc), 
so the expected increase in the number of edges per unit 
time is c — r(k) . At time t the expected number of vertices 
and edges are n = (l — r)t + N and m = (c — r{k))t + Mo, 
so in the limit t — ¥ oo the mean degree obeys 



(k) = 



2m 



2(c-r(fc)) 
1 - r 



and solving for (k) gives 



(k) = 



2c 



1 + r 

The average attractiveness is then 

oo 

(A) = A kPk = (k) + k* 

fe=0 



2c 



1 



(1) 



(2) 



3.2 Master Equation 

Let n be the number of vertices at time t. The expected 
number of vertices with degree k is — np^. One time 
step later this is n' k = (n + 1 — r)p' k where p' k is the new 
value of pk ■ Thus 

(n + 1 - r)p k = np k + 5 kc + — {A k _ 1 p k _ 1 - A k p k ) 

+ r(k + l)pk+i - rkp k - rp k . (3) 

The term Skc corresponds to adding of a vertex of degree c 
to the network. The term cA k _-yVk-\l (A) is the probabil- 
ity that a vertex of degree k — 1 gains an extra edge from 
the new vertex and becomes of degree k, and similarly 
cA k _ 1 p k ^i/ (A) is the flow from degree k to degree k + 1. 
The terms r{k + l)pk+\ and rkpk are the flows from k + 1 
to k and from k to k — 1 respectively, as vertices lose edges 
when one of their neighbors is removed from the network. 
Finally, rpk is the probability that a vertex of degree k is 
removed. Contributions from processes in which a vertex 
gains or loses two or more edges in a single unit of time 
vanish in the limit of large n and have been neglected. 

We are interested in the asymptotic form of the degree 
distribution pk in the limit of large t. Setting p' k = pk 
in (3) gives 



h,c + {Ak-iPk-i - A k p k ) 

+ r(k + l)p k+ i - rkp k - p k = . (4) 

as previously appeared in [17]. Naively, this equation ap- 
pears linear in the pk- But since it involves the mean at- 
tractiveness (A), the combination of (2) and (4) gives a 
nonlinear system of equations. 

When the attractiveness is proportional to the degree, 
Ak — k, the system (2) and (4) separates, and has been 
solved analytically in [17]. The authors showed that in 
this case the degree distribution exhibits a power-law tail 
in the case < r < 1 of growing networks, and follows 
a stretched exponential in the constant-size case r = 1. 
For more general attractiveness functions like the one in 
this paper, a fully analytic solution seems more difficult. 
Thus we focus on the behavior of pk for large k. Depending 
on the model parameters we find either a degree distribu- 
tion with a power-law tail (Figure 2) or an exponential 
tail (Figure 3). We confirm our calculations with numeri- 
cal simulation of the master equation, and through direct 
simulation of the network dynamics. 



3.3 High-Degree Expansion 

We now specialize the master equation to our model. Sub- 
stituting (2) into (4) gives, for k ^ {c — 1, c, c + 1}, 



In the case r — 1 where the network has constant size, we 
have (k) = c and (A) — c+ k*. 



(A) 



((*-l + fc*)p fc _x - (k + k*) Pk ) 



r{k + l)pk+i - rkp k -p k = 0. (5) 
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We will determine the asymptotic behavior of pk with 
a "high-degree expansion" , by approximating the ratio 
Pk/Pk—i as a Taylor series in 1/fc. We find a phase tran- 
sition between power-law and exponential behavior, and 
determine the phase diagram as a function of the param- 
eters c, r, and fc*. 

As an ansatz, assume that pk is a power-law times an 
exponential: 

p k = Ck a p k . (6) 

We can determine a and /? by taking fc 3> 1, and expand- 
ing the ratio Pk/Pk-i to leading orders in 1/fc. This gives 



Pfc-i V k 



O 



Pk+i 
Pk-i 



r i + 



2a 

T 



Substituting this into (5), multiplying by (A}/pk-i 
ignoring 0(1/ k) terms yields the equation 



(7) 
(8) 
and 



k(p-l)(r(A)p-c)- 



a(r(A)(3(2f3-l) 
-l) + c(k*(l-p) 



cp) 

1)=0. (9) 



Since (9) must be true for all fc, we can set the coefficient 
of fc to zero. This gives two solutions for /3, namely (3=1 
and 



P = 



-(A) 



(10) 



If r(A) > c, the solution p < 1 of (10) is physically 
relevant andp^ decays exponentially. However, if r(A) < c 
then (10) would give (3 > 1, which does not correspond to 
a normalizablc probability distribution. In that case (3=1 
is the relevant solution, and pk ~ k a is a power-law. Thus 
a phase transition occurs at r(A) = c. Applying (2), we 
can write this in terms of a critical value of fc* , 



fc: = 



c(l-r) 



(11) 



r(l + r) ' 

We illustrate the resulting phase diagram in Figure 1. 
3.4 The Power-Law 

To solve for the power-law exponent a, we again use (9), 
but now set the constant term (with respect to fc) to zero. 
If P = 1, this gives pk ~ fc" where 



(A)(l-r) + c 
c - (A)r 



< 0. 



(12) 



Note that a approaches — oo as we approach the transi- 
tion. As we show in Section 3.6, at criticality pk takes a 
stretched-exponential form. 

Substituting (1) and (2) into (12), we can express a in 
terms of c, r, and fc* as 



Q 



c(3-r) + fc*(l-r 2 ) 
c(l — r) — fc*r(l + r) 
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Fig. 1. Phase diagram of the degree distribution with 
c = 20. The two black dots indicate locations in the phase 
space for data shown on Figures 2, 3 and 4. 



In the special case fc* = 0, this recovers the result of [17] 

3-r 
1 — r 

while in the special case r — 0, it includes the result of [10] 

fc* 

a = — 3 . 

c 

Finally, since fc* > and r € [0, 1], we have 

— co < a < — 3 , 

so that the degree distribution has a finite average as well 
as (except when fc* = r = 0) a finite variance. 

Above the transition, ft is given by (10). Again setting 
the constant term of (9) to zero gives 

C (fc*-l) + (A)(l + r(l-fc*)) 

a = TT\ • 13 ) 

[A)r — c 

Thus we have a power-law correction to the exponential 
decay, pk ~ fc a /3 fe . In terms of our parameters, 



a = fc* — 



c(l + r) 



r(2c + fc*(l + r)) 

(l + r)(c + fc*(l + r)) 



fc*r(l + r) - c(l - r) 
Note that in this regime a may be positive. 

3.5 Finite- Degree Corrections 

We have derived the leading behavior of pk for large fc, 
namely a power-law times an exponential. In this section, 
we obtain the next-order correction, by taking the Tay- 
lor series to second order in 1/fc. This correction becomes 
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important in the exponential regime, where the exponen- 
tial decay of the degree distribution makes smaller degrees 
more relevant. 

Rather than starting with an ansatz for the correction 
term, we derive it by expanding Pk/pu—x to second order 
in 1/k. Write 



(14) 




! ^alnk + {a + a l /2 - K +0 (l)) . ,1.5, 

where ~ hides multiplicative constants. Here we used the 
Taylor series for the logarithm, 

ln(l + e) = e-e 2 /2 + 0(e 3 ), 

the approximation for the fcth harmonic number 

k 

£i=lnfc + 7+ ^+0(l/fc 2 ), 

2 — 1 

where 7 is Euler's constant, and 



E 



1 IT 2 1 



- + 0(l/fc 2 ) 



Eq. (15) gives a multiplicative correction of the form 
c s / k to our earlier form (6) for p^, 



where 



Pk = C(i k k a e 5 ' k (1 + 0{l/k 2 )) 



- a + a 

= K . 

2 



To determine k and therefore 5, using (14) we write 



Pk+i 
Pk-i 



P 1 + -r + 



2a 2k + a 2 



k 2 



O 



(16) 
(17) 

(18) 



Substituting (14) and (18) into the master equation (5) 
and multiplying by (A)/pk-i, we obtain an equation akin 
to (9), with terms of order fc, 1, 1/k, and negligible terms of 
order 0(l/fc 2 ). The coefficient of k and the constant term 
are identical to those in (9), giving the same solutions for 



a and /3 as before. Setting the coefficient of 1/k to zero 
and applying (17) gives 



6/3 - 4a/3)) 



a c(l + a + k*) + (A) (2 + r(l + a 
~ 2 c + r(A)(l-2P) 



Multiplicatively speaking, the correction term e s / k be- 
comes negligible as k — > oo. However, it makes a signifi- 
cant difference for small values of k, and greatly improves 
agreement with the simulations in the next section. We 
note that the same technique, expanding the ratio between 
Pk-i, Pk, and Pk+i to higher degree in 1/k, can give us as 
many correction terms as we wish. 



3.6 The Stretched Exponential At Criticality 

We saw above that as we approach the critical point, the 
power-law exponent a diverges to — oo, and the exponen- 
tial factor p approaches 1. In this section we show that 
the degree distribution in fact becomes a stretched expo- 
nential at this point, due to the appearance of half- integer 
powers of 1 /k in the high-degree expansion of Pk/Pk-i- 
We start with the ansatz 



Pk 



Ck a fi k e cVk . 



Expanding to order k 3 / 2 , we have 



Pk 
Pk-i 
Pk+i 
Pk-i 



= 01 



c a + e 



= /3 2 1 



2-sfk 

v Vk 



k 

2a + C/2 



C 3 /48 + aC/2 + C/8 

fc 3/2 

C 3 /6 + 2c< 

fc 3/2 



with error terms of order 1/fc 2 . Substituting this into the 
master equation (5) as before, if ( ^ then the terms of 
order k and Vk force b — 1 and c = r(A). In other words, 
£ can be nonzero only at the critical point. The term of 
order 1 then gives 

C = -2/Vr, 

and the term of order /c -1 / 2 gives the power-law correction 

3 k* 
a = 1 . 

4 2 

This recovers the results of [17] for the special case 
fc* = and r — 1, where £ = —2 and a — —3/4. However, 
these calculations are significantly more technical, evalu- 
ating generating functions and their derivatives in terms 
of special functions. Our high-degree expansion is closer in 
spirit to [13], where pk is written as a telescoping product 
of ratios Pk/pk-i- 



4 Simulations 

To check that our asymptotic calculations are correct, we 
conducted two kinds of simulations: direct simulation of 
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Fig. 2. Comparison between simulations for k* = 10, 
r = 0.10, c = 20 and our solution, which gives a power- 
law pk ~ k a with a = —4.02. There is good agreement for 
k > 200. This set of parameters is in the power-law regime, 
below the phase transition at k* 163.6. Note the log- 
log scale. The integrated master equation dips down at 
k « 3000 because it has not yet reached equilibrium at 
the highest degrees. 

We show a network simulation for a single network. For 
large degrees, we use logarithmic binning, so that each 
point represents the average over an interval of degrees of 
width proportional to logfc. 



the dynamics of finite networks, and numerical integration 
of the master equation. 

In our direct simulations, we grew the network stochas- 
tically according to the model, up to size n — 5 x 10 7 . How- 
ever, it is hard to explore the tail of the degree distribution 
in the exponential regime, since pk falls off exponentially. 
For instance, to measure a probability p^ ~ 10~ 10 we 
would need a network of size more than n = 10 10 , unless 
we use large bin sizes. In this case, numerically integrating 
the master equation (3) until it reaches equilibrium lets us 
explore the asymptotics of pk far more efficiently. 

Since the normalization constant C depends on the 
values of pk for small fc, which our analysis does not try 
to predict, we adjust C to fit the simulations. We do not 
tune any other parameters. In particular, a and /3 are 
determined by our analysis, rather than fit to the data. 

Figure 2 shows results in the power-law regime below 
the transition. There is good agreement between our so- 
lution and both types of simulations above k > 200 or so. 
The direct simulation differs somewhat from the master 
equation at large k due to finite-size effects. 

Figure 3 shows results in the exponential regime, above 
the phase transition. Here as well, there is good agreement 
between simulations and our asymptotic solution for large 
enough k. Figure 4 shows the same parameters at larger 
degrees; as discussed above, we reach these larger degrees 
by abandoning direct simulation and integrating the mas- 
ter equation. The agreement with our asymptotic solution 
is excellent. 
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Fig. 3. Comparison between simulations for k* = 10, r = 
0.80, c = 20 and our solution, which gives pk ~ ft k k a e s / k 
with /3 = 0.776, a = 3.423 and S = 80.89. There is good 
agreement for k > 30. This set of parameters is in the 
exponential regime, above the phase transition at k* ~ 
2.78. Note the semi-logarithmic scale. 
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Fig. 4. Comparison of our analytic solution and the in- 
tegrated master equation for k* = 10, r = 0.80, c = 20 
(the same parameters as Figure 3) at larger degrees than 
direct simulations can reach. 



The pseudo-code and source code for the simula- 
tions can be found at http://www.rouquier.org/jb/ 
research/papers/2010_growing_network/. 



5 Conclusion 

We have studied dynamical networks that are generated 
by a model where growth takes place through a combi- 
nation of preferential and uniform attachment, and where 
nodes are removed randomly at a certain rate. Both growth 
and node removal are key features of many real- world net- 
works. Nodes in peer-to-peer networks may be added or 
removed, people join and leave social networks, nodes in a 
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communication network can be attacked or degrade with 
time, and so on. Uniformly random attachment appears, 
for instance, if people choose random seats and strike up 
conversations with their neighbors, or are assigned to ran- 
dom classrooms; random attachment also occurs, by de- 
sign, in some peer-to-peer protocols. 

We have solved for the asymptotic degree distribution, 
and found a phase transition between power-law and ex- 
ponential behavior, with a stretched exponential at the 
critical point. Thus, in contrast to pure growth models, 
an asymptotically linear attractiveness function is not a 
sufficient condition for a power-law degree distribution. If 
the growth rate is too small, or the node removal rate is 
too high, the degree distribution is exponential. 

Our findings are relevant for real-world networks where 
both growth and node removal are important. They also 
imply further potentially interesting consequences for the 
evolution of networks whose effective growth rates vary 
with time; for instance, for networks that start with a 
high growth rate, but that reach a state where nodes are 
added and removed at about the same rate, e.g. due to 
limits on the network's overall size or population. 

Similarly, it would be interesting to understand how 
the macro-dynamics of a network change as the growth 
and/or removal rates are varied so that we approach or 
cross the topological phase transition. This includes the 
approach to the asymptotic degree distribution from a 
noncquilibrium initial state; whether this approach shows 
critical slowing down near the transition; and the dynam- 
ics after a sudden change of the growth and/or removal 
rates, say from a region from the power-law regime to the 
exponential one. 

A broader question is how the transition affects various 
types of dynamics taking place on the network, such as 
search [18], congestion, and robustness to attack. 

Finally, another direction for future work is to intro- 
duce some kind of quenched disorder into the network 
model. This could include allowing fc* to vary from node 
to node, based on the node's intrinsic "fitness" or "attrac- 
tiveness" [14,19]. We believe that the asymptotic behavior 
of the degree distribution and its phase diagram will be 
similar to our results here as long as k* has bounded ex- 
pectation and variance. 
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